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Turbulent Flow in Shock Tubes of Varying 

Cross Section* 

Robert F. Dressier 

The unsteady escape flow of a compressible gas is investigated subject to the influences 
of varying duct cross section and mechanical retardation due to turbulence and frictional 
dissipation. This study is applicable to shock tubes, intermittent wind tunnels, and general 
pneumatic-control systems, such as airbrakes. Wave-front velocities and trajectories are 
obtained by direct integration along the bounding Mach lines. In the general case, calcu- 
lation of first-order approximations is reduced to a set of independent ordinary differential 
equations for easy stepwise integration. For the special cases of constant resistance or 
exponentially converging-diverging ducts, first-order solutions are obtained analytically. 
An error estimate is derived. Local Mach number variations are presented throughout the 
flow. 

1. Introduction 

One type of wind tunnel in present use employs the principle of intermittent discharge of 
atmospheric air into a vacuum chamber, or intermittent high-pressure discharge into an 
atmospheric chamber. In the latter case, a rapidly advancing shock discontinuity will also 
he generated in the forward wave. The familiar Kiemann shock tube furnishes such an 
example; for this type of duct Riemann [l] 1 studied the basic theory of one-dimension a I unsteady 
gas flows. In the present paper the emphasis will be upon the flow immediately after the 
diaphragm is broken, not upon the flow developing after a sufficiently long time when an almost 
steady-state regime is approached. The intermittent-vacuum-type system gives a more con- 
stant stagnation temperature and requires less power than the intermittent-pressure type, but 
gives lower densities in the test section. In the first stages of operation, these flows in shock 
tubes or vacuum-type wind tunnels can be approximately described by the known centered 
simple-wave solution [2, ch. Ill], defining a subsonic-supersonic pattern about the diaphragm 
position. It is the purpose of this paper to investigate quantitatively the possible improvement 
in density, velocity, or Mach number that might be created if such an unsteady flow is made 
to pass through a test section of varying cross section; at the same time the retarding effect 
due to wall resistance and ensuing turbulence will also be considered. 

Any forward shock discontinuity that would be created in an actual flow if it propagates 
into a region containing some gas will be ignored in the present treatment. Even if a relatively 
strong shock exists, the results to be derived here will still be approximately applicable to the 
backward moving rarefaction wave region. This will furnish quantitative information, for 
example, concerning frictional effects in the rarefaction waves propagating through pneumatic- 
control circuits, such as airbrake lines. 

The case of steady flows in variable ducts has recently been studied by Meyer [3], empha- 
sizing the propagation of advancing and receding wave fronts due to a finite disturbance super- 
imposed upon a basic steady flow, with consideration given to the varying cross-section effect, 
but not the resistance effect. 

Dimensional quantities (barred) will first be introduced to define the dimensionless ones 
(unbarred), which will then be used throughout this paper. In the standard notation as given 
in [2], x is taken in the direction of flow; time is denoted by t, with the valve or diaphragm 
at x = opened at £=0. Let p denote pressure, p the mass density, and u the velocity com- 
ponent in the downstream direct ion. We consider the nonlinear flow of a poly tropic gas having 
the equation of state 
_____ P=Ap\ (1) 
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Figure 1. 



with known constant A and adiabatic exponent 7. We introduce in the usual manner the local 
sound speed c into the differentiated equation of state, 



c 2 =lA~p y ~ l . 



(2) 



For F<0, the gas is at rest to the left of the valve, x<0, and has properties p , c , andi2 =0. 
For F<0, there is no gas in x>0 (see fig. 1). We assume the flow to occur in a duct with cross 
section varying gradually enough to keep the flow essentially one-dimensional with uniform 
velocity over any cross section, and subject to the retarding force from the wall friction 
and turbulence, depending empirically upon the square of the velocity and the roughness 
coefficient for the wall surface. The dynamic equation is then 



—+u ^=H c ■^=+eR(x)u 2 =0, 

&t da: 7—1 bx l 



(3) 



where e is a small dimensionless expansion parameter, and the given function R(x)^>0 varies if 
the roughness of the wall surface changes from region to region, or if the cross section varies, 
thus varying the hydraulic radius 2 of the duct in a known manner. Letting 3(z) denote the 



area of the cross section of the duct, the modified equation of continuity will be 



bt ^ 



«,... bu , , 



■l)ea(x)uc=0, 



(4) 



in which a(x) = [doo/dx]/oo> with ea = corresponding to a uniform duct. 

In eq (1), (3), and (4) we have tacitly assumed constant entropy for the flow. This is a 
simplification tantamount to neglecting the heat change (and hence its extra resulting effects 
upon the mechanics of the gas flow) which will always be produced by the operation of the 
resistance term in converting kinetic energy into thermal energy. The present discussion will 
consider the mechanics ofj the flow including only [the retarding force, not the full thermo- 
dynamical aspects of the problem. To do so would require also the inclusion of the energy 
equation, creating a system of three equations in the three unknowns: particle velocity, sound 
velocity, and specific entropy. Such a system has an associated set of six characteristic equa- 
tions; these have been used by Hall [4] and others for purposes of quite arduous numerical 
integrations performed over characteristic meshes in the (x,t) plane. The simplified model 
used in our present paper can be described also in the following manner. We consider a one- 
dimensional flow through a tube with the mass traveling [in thin slices perpendicular to the 
walls. Each slice is subjected to a retarding shear force distributed over its peripheral area. 
This force F(p,u) = —\pu 2 varies with the changing density and velocity of the moving slice. 
Results obtained for this model will thus extract the purely mechanical effect of the friction 
from the full nonisentropic thermodynamical effects. 



2 If k(x) measures the roughness of the wall material, then R= k(x)fr(x), where the hydraulic radius r is defined as w/>. P is the perimeter of 
the cross section having area 
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Now introduce the unbarred dimensionless quantities defined by x=(g/cl)x, t = (gfc )t, 
u=(l/c )u, c=(l/c )c,p=(l/p )p, R=(cl/g)R, a=(c 2 Jg)a, and abbreviate 2/(7—1)=/, the un- 
modified escape speed. Equations (3) and (4) then become 



bu . bu . 7 be . D / N 2 n 
Zi +U b-x +lc Zx +tRix)u=0 

7 be . 7 be bu t N A 
Z ^-.+m =- — \-c ^ — [-€a{x)uc=0. 
bt bx bx w 

The escape flow solution for €=0 (uniform duct and resistance neglected) is 

2 



u v 



~7+l 

.7—1 
"7+1 



(1+m) 



(Z-m), 



(5) 
(6) 

(7) 
(8) 



where m=x/t. This is supersonic in the forward region 0<Cx<lt; subsonic in the backwave 
— t<x<i0; and u°=e°= 2/(7+1) at the fixed sonic point x=0. 
The equations for the Mach lines of the system (5), (6) are 

T dx 



T dx , 



II + :d(u—lc) = —e(Ru—ac)udt 
II-\d{u-\-lc) = —e(Ru-\-ae)udt, 



(9) 



where the /+ and II + relations hold along the b+ family of receding Mach lines in the (x,i) 
plane, and the other pair along the £_ advancing lines. The Mach lines OB and OF in figure 
2 define the moving limits of the flow zone. As the quantities R{x) and a(x) do not appear in 
the first equation of (9), it follows that the curve OB must remain a straight line with direction 
m= — l in all cases. 
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2. Locus of the Wavefront 

It is seen at once from the last two equations in (9) that, within the accuracy of the basic 
relations (5) and (6), the duct shape has no effect upon the wavefront trajectory where c = 
because the terms containing a(x) will vanish along this wave front. Its locus is a function of 
eR, however, and will remain a double Mach line as the plus and minus relations coincide for 
c = in (9). The / and II pair of equations can then be integrated explicitly to obtain the 
exact wavefront locus, which is 

tFix) =^j\f>™ dt , (10 > 

in which e is included in the notation eR=dt. The velocity of the wavefront is given by 

M*)-^™*. (ii) 

In particular, if 9J = constant, e. g., in a duct having homogeneous wall surface and uniform 
cross section, the position and velocity will be 

(12) 



These exact results will be used later to furnish an error estimate for the approximations to 
be derived for the interior wave zone. The forward wave velocity will not approach asymp- 
totically any nonzero value, but will continually slow down, approaching, but never attaining, a 
state of rest. This is true more generally for any variable R(x) with a(x) arbitrary, as can be 
seen from (11), because R(x) is always nonnegative. 

The Mach lines inside the flow will be locally distorted by the influence of R and a from 
their positions for the basic solution (7) and (8), which are 

b+: x=mt, -l<m<2/(7— 1), (13) 

where each 6_ curve emanates from the point ( — p,p) on OB. 

3. Comparison with Hydrodynamic Problem 

If e = in (3) and (4), the equations become identical for 7 = 2 with the approximate 
Saint-Venant equations for two-dimensional water flow with a free surface, where c 2 is the 
water height, and u is the horizontal velocity of a cross section. If this analogy extended 
also to the equations when the resistance term is included, then the present aerodynamic 
problem would be analogous to the hydraulic "dam-break" problem previously considered by 
the author in [5]. The hydraulic-resistance function, however, is more complicated in an 
essential way, taking the form — Ru 2 /c 2 rather than — Ru 2 . The hydraulic dam-break wavefront 
contains a second-order pole singularity, and its locus is thereby not a characteristic curve 
itself, but only an envelope of such curves. Explicit integration for its trajectory is apparently 
not possible in that case. The simpler nature of the present problem permits carrying the 
different mathematical results somewhat farther and with less approximation involved for the 
basic model used. 



256 



4. Functions for First-Order Effects 

Using- (7) and (8), we substitute 

u(z f t,€) = u°(x,t) + U(x,t)e + 0(e 2 ) 

(15) 
c(x,t,e) =c°(x,t) + C(x,t)e + 0(e 2 ) 

into (5) and (6). The terms in e define the system 

< 7 + ] ) al + 2 (7 +1 ) ai +2 (/t) 5,+i ^"i C+ t+i (j +1 ) B( * )==0 

(16) 

(7 _„( ; _ f )^ +;( , +0 ^ +a(f+1) g_ ?[;+?c+2( ^ )(( _ f)(f+I)(lw=0 

This system for [7 and C has the same known Mach lines in the (x,t) plane as defined in 
(13) and (14), along which the following Riemann relations will hold: 

b + : d(U-W)=\j-^ t {C-U)-{Ru Q -ac")u^\ di 

(17) 
b_: d(U+lC) = -(Ru +ac )v dt. 

This second equation is independent of U and C on the right side; therefore, it is possible 
for given F(x) and a(x) to integrate the line integral explicitly along each known 6_ curve, 
using (14). The result will be a known expression of form 

U+lC=$(t } p). 

This will then permit a very easy separate step-wise numerical integration along each 6 + 
curve to solve for U and C. Thus the problem is reduced to a set of independent "one-dimen- 
sional" numerical calculations for a single ordinary differential equation along each b + curve. 
This is a calculation of much greater simplicity than the two-dimensional numerical procedure 
that would be necessary for the full problem in (5), (6), and (9). 

5. Integration in Special Cases 

Instead of using the numerical approach for arbitrary R and a, we return to solve the basic 
eq (16) explicitly for the special case when R and a are both constants. In particular, a = 0, 
R= constant >0 will describe the resistive effects on the flow in a duct with an3 r fixed cross 
section and homogeneous wall roughness. The combination a = constant 5^0, R=0 will describe 
the effect of a converging or diverging duct, neglecting resistance effects. For example, a>0 
can describe a duct with circular cross section, exponentially diverging, or with rectangular 
cross section, constant width, also exponentially diverging; a<^0 corresponds to exponentially 
converging ducts of either type. The combination R and a, both nonzero constants, is less 
meaningful physically, because in such a varying duct the wall roughness would have to be 
varied in conformity with the changing hydraulic radius in order to maintain R constant. 
Therefore, the contributions of R and a will be considered separately, although the solution 
will be made for both simultaneously. 

First, we put 7=1.4 for air to shorten the expressions; the escape velocity is then l=b. 



257 



System (16) has for its general solution 
C7fe , ) =_^0)_-|[ 8 4_36(5-f)+5(5-f) 2 ]i?^| § [2(5-|)-(5-|) 2 ]a i +ir( a ) J 

C(x,0=^(|)- I | i [60-36( 5 -f)+7(5-|) 2 ]^-^[lo( 5 -f)-3(5-f) 2 ]^+|+ (18) 

where a={b—xjt)t llz , K and F are arbitrary functions, and 5 is an arbitrary constant. The 
form of this solution indicates that U and C must be either infinite or zero at 2=0; hence we 
must choose K=0 and 5=0 to keep the functions zero. The boundary conditions U=C=0 
along line OB then determine the function F=(l/1296) (50R— 15a) (5—x/t)H. With x/t=m, the 
main results reduce to 



or: c 

U=~-^- (2m 3 +l5m 2 +24m + ll)Rt+-^- (m 3 -6m 2 +3m + 10)<^ 
C=~—^— (4m 3 +3ra 2 -6m-5)i^+ -^ (2m 3 -3m 2 -30m-25)a*. 



(19) 



Denoting eRt=<r and eat=r, the flow velocity and local sound speed correct to first order can 
be written from (19), using (7) and (8), in the form 



u=u° -\-hi(m)<T-\-h 2 (m) r 
c=c°+ki{m)<r J rk2(rri)T. 



(20) 



Once u and c have been calculated by (19) and (20), the pressure and density are immediately 
known through relations (1) and (2). 

5.1. Effects on Velocity and Density 

To provide immediate perception of the various effects, the h, k polynomials are graphed in 
figure 3, in which hi is plotted to a separate vertical scale given on the right. We note that an 
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exponential duct causes an antisymmetrical effect on the velocity about the moving midpoint 
of the flow zone (m=2), whereas the effect on sound speed and density is almost symmetrical 
about m = 2.8. On the other hand, in a uniform duct, the resistance effect is a monotonic 
decrease in the velocity as it approaches the wavefront ; it creates a very slight increase in sound 
speed at the rear of the flow, then a larger decrease in the forward region. 

As a numerical example to illustrate the various separate effects, u, c, and the local Mach 
number M=u/c are evaluated in figure 4, A, for o-=l/10; and for r=+l in 4, B, and r= — 1 
in 4, C. The dotted lines represent the quantities u°, c°, and M° when a = r = 0. 

5.2. An Error Estimate 

As shown by the c curve in figure 4, A, the abscissa of the point where the solid curve inter- 
sects the m axis defines the position of the wavefront to the first order approximation. Like- 
wise the corresponding value of u on the curve directly below furnishes an approximation for 
the velocity of the wavefront. We have, however, already obtained the exact expressions 
for the position and velocity of the wavefront in (12). The approximations and the full solu- 
tion agree exactly along line OB; the resistive effect is a cumulative one with increasing velocity 
in the forward part of the wave; hence it is to be expected that the maximum of the approxima- 
tion- error will occur at the wavefront. Thus the error bound obtainable in this way should 
also be valid for all portions of the wave. Such a calculation for the position of the wavefront 
indicates that the approximate result will be about 10 percent too large when c grows as large 
as 0.2. The other case of variable cross section does not permit this simple method for an 
error estimate. 

5.3. The Sonic Path 

If x=s(t,e) is the path along which u=c, with s=0+#€+0(€ 2 ), the resulting double series 
expansion correct to first order gives the relation S(t) = [C(0,t) — U(0,t)]t. For the case R= 
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Converging Duct 



constant, a = 0, this defines the sonic path correct to first order to be the downstream parabola 

x=0.23eRt 2 . 

For the other case a = constant, R=0, we get the parabola 

.7 = — 0.18ea£ 2 , 

going upstream or downstream as the duct is diverging or converging, respectively. 

5.4. Flow Past a Fixed Point 

The Mach number distribution across the entire flow at a fixed time has been presented in 
figure 4. Of more practical interest, however, from the standpoint of applications in shock 
tubes or wind tunnels is the time variation of local Mach number at a fixed point downstream 
in the duct, at a point where a fixed model might be placed. This can be derived in an obvious 
manner from the expressions already presented. As an illustrative example, the Mach number 
is shown in figure 5 at the fixed point £=100 for eR=0.005, t = in 5, A; for ea = +0.01, or = 
in 5, B; and for ea=— 0.01,(r = in 5, C. The dotted curve in each figure again represents the 
M° corresponding to the undisturbed solution. 
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Figure 5. 
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